#include "dyn_specialFunctions_besseli.h"

BesselI::BesselI(double _order)
{
    setOrder(_order);
}

double BesselI::operator ()(const double x)
{
    if (order < 0.0)
    {
        double Pi = 3.1415926535897932384626433832795028841971693993751058;
        return boost::math::cyl_bessel_i(-order,x)+(2.0/Pi)*sin(-Pi*order)*boost::math::cyl_bessel_k(-order,x);
    }
    return boost::math::cyl_bessel_i(order,x);
}

d_complex BesselI::operator ()(const d_complex z)
{
    // Input values
    double zr = real(z);
    double zi = imag(z);
    double nu = fabs(order);
    int kode = 1;
    int N = 1;

    // Output values
    double cyr,cyi;
    int nz,ierr;

    zbesi_(&zr,&zi,&nu,&kode,&N,&cyr,&cyi,&nz,&ierr);

    d_complex answer(cyr,cyi);

    if (order < 0.0)
    {
        double cyr2,cyi2;
        double Pi = 3.1415926535897932384626433832795028841971693993751058;

        zbesk_(&zr,&zi,&nu,&kode,&N,&cyr2,&cyi2,&nz,&ierr);

        d_complex answer2(cyr2,cyi2);

        return answer+(2.0/Pi)*sin(Pi*nu)*answer2;
    }
    return answer;
}
